## Backup for statements in draft
# Emissions reductions, US vs global
p.base = sim.base$sim$global.sim[year(date)>=2020]
p.ban = sim.ban.lag.0$sim$global.sim[year(date)>=2020]
p.18 = sim.18.75.royalty.rate.on$sim$global.sim[year(date)>=2020]
deltas.sep$sim.18.75.royalty.rate.on$deltas[.N, (US.Oil-Fed.Oil)]
deltas.sep$sim.18.75.royalty.rate.on$deltas[.N, (US.Oil-Fed.Oil)/-Fed.Oil]
deltas.sep$sim.18.75.royalty.rate.on$deltas[.N, (US.Gas-Fed.Gas)/-Fed.Gas]
deltas.sep$sim.18.75.royalty.rate.on$deltas[.N, (US.Oil-Fed.Oil)/-Fed.Oil]
deltas.sep$sim.18.75.royalty.rate.on$deltas[.N, ((US.Oil-Fed.Oil)*6+US.Gas-Fed.Gas)/-(Fed.Oil*6+Fed.Gas)]

# BAU average oil and gas royalties are projected to be $9-10 billion/year
gov.revenue(sim.base)$roy.rev.summary[year(date)>=2020, mean(total.revenue)]
gov.revenue(sim.ban.lag.0)$roy.rev.summary[year(date)>=2020, mean(total.revenue)]
deltas$sim.ban.lag.0$deltas[.N, Global.Emis.Tot]

# As a point of comparison, the emissions associated with federal (total) US oil and gas 
# production amounted to about 719 (4,622) million tons of CO$_2$e in 2019.
# EIA: 2019 oil and gas production was 12.248 mb/d and 112.034 bcf/d (40,892,458 MMcf)
# ONRR: 2019 oil and gas production was 2.718665 mb/d (697651763+294660842)/365/1e6 
# and 12.11958 bcf/d (1061196836+3362450982)/365/1e6
# This means that federal & total O&G emissions were
2.718665*365*oil.emis.rate.fed + 12.11958*365*gas.emis.rate.fed
12.248*365*oil.emis.rate.nonfed + 112.034*365*gas.emis.rate.nonfed
# https://www.eia.gov/dnav/pet/pet_crd_crpdn_adc_mbblpd_a.htm
# https://www.eia.gov/dnav/ng/ng_prod_sum_a_EPG0_FGW_mmcf_a.htm
# https://revenuedata.doi.gov/?tab=tab-production

# In other words, about 30 percent of the leakage arises from US production 
# from nonfederal lands and the other 70 percent from foreign (ROW) supply.
sum.table.num[,'Nonfed US emis']/(sum.table.num[,'Nonfed US emis']+sum.table.num[,'ROW emis'])

# under a leasing ban, "the emissions reductions from crude oil are about twice as large as those from gas."
deltas.sep$sim.ban.lag.0$deltas[.N, Global.Emis.Oil/Global.Emis.Gas]
deltas.sep$sim.carbon.adder.fed.20$deltas[.N, Global.Emis.Oil/Global.Emis.Gas]

# estimated climate benefits of approximately \$5-\$8 billion annually for 
# a carbon adder and \$6-\$10 billion annually for a leasing ban. 
load(working.data.loc)
# load(working.data.loc.high.e)
carb.add.compare = calc.deltas(sim.base, sim.carbon.adder.fed.20, yrs=2020:2050)
carb.add.compare$deltas[-.N, year.num:=as.numeric(year)]
carb.add.compare$deltas[-.N, scc := 50*(1+0.02)^(year.num-2020)]
carb.add.compare$deltas[-.N, bene := scc*-Global.Emis.Tot/1e3]
carb.add.compare$deltas[-.N, mean(bene)]

# Global versus US Gas markets
# Only after 2045 or so do exports begin to exceed 25 bcf/d, surpassing capacity that already exists or is under construction.
dev.off()
exports = sim.base$sim$global.sim[, .(date, US.gas.supply, US.gas.demand)]
exports[, gas.export := US.gas.supply - US.gas.demand]
exports[, mean(gas.export), by=year(date)]

load(working.data.dir%&%'/working_fed_lands_data_modelrun_Domestic_gas_market_2021-11-01.RData')
sum.table.num
load(working.data.loc) # reload base case data

deltas$sim.carbon.adder.fed
deltas.summary$sim.carbon.adder.fed

tas(sim.base, sim.ban.lag.0, yrs=c(2020:2030))$deltas[, .(year,Fed.Emis.Tot, Global.Emis.Tot, 1- Global.Emis.Tot/Fed.Emis.Tot)]
calc.deltas(sim.base, sim.ban.lag.0, yrs=c(2020:2030))$deltas[-.N, 1- sum(Global.Emis.Tot)/sum(Fed.Emis.Tot)]

# Because US federal oil and gas production is a relatively small fraction of global
#  oil and gas supply (about 2 percent for oil and 3 percent for gas)
onrr.prod[year(date)==2018, mean(Oil_Federal_Offshore+Oil_Federal_Onshore)]/
  weo.projections[year(date)==2018, mean(world_oil_cons)]

onrr.prod[year(date)==2018, mean(Gas_Federal_Offshore+Gas_Federal_Onshore)]/
  weo.projections[year(date)==2018, mean(world_gas_cons_bcfd)]

# Federal oil and gas production as a % of global
sim.base$sim$us.sim$OilProduction_AllWells[year(date)==2018, mean(Oil_Federal_Offshore+Oil_Federal_Onshore+Gas_Federal_Offshore+Gas_Federal_Onshore)]/
  sim.base$sim$global.sim[year(date)==2018, mean(Tot.Oil.Supply)]
sim.base$sim$us.sim$GasProduction_AllWells[, mean(Oil_Federal_Offshore+Oil_Federal_Onshore+Gas_Federal_Offshore+Gas_Federal_Onshore)]/
  sim.base$sim$global.sim[year(date)==2018, mean(Tot.Gas.Supply)]

# "Leakage to US producers on state and private lands constitutes about one-third of the total 
# leakage effect,
deltas$sim.ban.lag.0$deltas[year=='Cum. Avg.', (US.Emis.Tot-Fed.Emis.Tot)]/ # US nonfederal leakage
  deltas$sim.ban.lag.0$deltas[year=='Cum. Avg.', (Global.Emis.Tot-US.Emis.Tot)] # non-US leakage
deltas$sim.carbon.adder.fed.20$deltas[year=='Cum. Avg.', (US.Emis.Tot-Fed.Emis.Tot)]/ # US nonfederal leakage
  deltas$sim.carbon.adder.fed.20$deltas[year=='Cum. Avg.', (Global.Emis.Tot-US.Emis.Tot)] # non-US leakage

# despite the fact that those sources historically represented less than 15
# percent of global oil and gas supply."
# Nonfederal oil and gas production as a % of global, in BOE
(sim.base$sim$us.sim$OilProduction_AllWells[year(date)==2018, mean(Oil_Nonfederal_Offshore+Oil_Nonfederal_Onshore+Gas_Nonfederal_Offshore+Gas_Nonfederal_Onshore)]+
    sim.base$sim$us.sim$GasProduction_AllWells[year(date)==2018, 1/6*mean(Oil_Nonfederal_Offshore+Oil_Nonfederal_Onshore+Gas_Nonfederal_Offshore+Gas_Nonfederal_Onshore)])/
  (sim.base$sim$global.sim[year(date)==2018, mean(Tot.Oil.Supply)]+
     sim.base$sim$global.sim[year(date)==2018, 1/6*mean(Tot.Gas.Supply)])

#  only about 18.5 percent of baseline federal emissions over the 30-year 
# simulation horizon are from newly drilled wells during this 10-year window. 
if (FALSE) {
  sb = simulate.policy() # run this with the following "# for testing" lines uncommented,
  (sum(res$OilProduction_NewWellsCovered[year(date)>=2020, c(4,5,8,9)]*365.25/12*oil.emis.rate.fed)+
      sum(res$GasProduction_NewWellsCovered[year(date)>=2020, c(4,5,8,9)]*365.25/12*gas.emis.rate.fed))/
    (sum(res$OilProduction_AllWells[year(date)>=2020, c(4,5,8,9)]*365.25/12*oil.emis.rate.fed)+
       sum(res$GasProduction_AllWells[year(date)>=2020, c(4,5,8,9)]*365.25/12*gas.emis.rate.fed))
  
  # Check:
  sum(sim.base$sim$us.sim$OilProduction_AllWells[year(date)>=2020, c(4,5,8,9)]*365.25/12*oil.emis.rate.fed)+
    sum(sim.base$sim$us.sim$GasProduction_AllWells[year(date)>=2020, c(4,5,8,9)]*365.25/12*gas.emis.rate.fed)
}

# These elasticities start at about +0.2 in the short run (2020) for both oil and gas supply and rise gradually in the long run 
# (2050) to +0.9 for oil and +1.2 for gas, 
# the ROW oil and gas supply elasticities are +0.4 for oil and +0.5 for gas over the simulation horizon
weo.projections[year %in% c(2020,2050), mean(row_crude_elasticity, na.rm=T), by=year]
weo.projections[year %in% c(2020,2050), mean(row_gas_elasticity, na.rm=T), by=year]
weo.projections[year %in% 2020:2050, mean(row_crude_elasticity, na.rm=T)]
weo.projections[year %in% 2020:2050, mean(row_gas_elasticity, na.rm=T)]

## Graph of emissions and revenue results
load(working.data.dir%&%'/working_fed_lands_data_modelrun_high_elasts_2021-11-01.RData')
sum.graph.high = mapply(FUN=extract.vars.for.table, d=deltas[sims.for.table], nm=sims.for.table)
sum.graph.high = t(sum.graph.high)
load(working.data.dir%&%'/working_fed_lands_data_modelrun_2021-11-01.RData')
sum.graph.base = mapply(FUN=extract.vars.for.table, d=deltas[sims.for.table], nm=sims.for.table)
sum.graph.base = t(sum.graph.base)

sum.graph.base = sum.graph.base[c(3,4,6),]
sum.graph.high = sum.graph.high[c(3,4,6),]

sum.graph.base = cbind(sum.graph.base, 'Nonfed emis' = rowSums(sum.graph.base[,c('Nonfed US emis','ROW emis')] ))
sum.graph.high = cbind(sum.graph.high, 'Nonfed emis' = rowSums(sum.graph.high[,c('Nonfed US emis','ROW emis')] ))
pol.nms = c('25%\nRoyalty\nRate', '$50\nCarbon\nAdder', 'Moratorium')

sum.graph = sum.graph.base
palette(c(rgb(136,196,244, maxColorValue=255), # RFF mid-blue
          rgb(116,100,94, maxColorValue=255), # RFF brown
          rgb(80,177,97, maxColorValue=255), # RFF green
          rgb(4, 39, 60, maxColorValue=255), # RFF black
          rgb(224, 228, 231, maxColorValue=255)))
my.cols


# For both the carbon adder and moratorium policies, the majority of the estimated global emissions reductions come from reduced oil consumption
deltas$sim.carbon.adder.fed$deltas[4,'Global.Emis.Oil']/deltas$sim.carbon.adder.fed$deltas[4,'Global.Emis.Tot']
deltas$sim.ban.lag.0$deltas[4,'Global.Emis.Oil']/deltas$sim.ban.lag.0$deltas[4,'Global.Emis.Tot']

###### Model validation plot
spud.sim = sim.base$sim$us.sim$spuds
spud.sim = spud.sim[,.(date,
                       Oil = Oil_Nonfederal_Onshore +  Oil_Federal_Onshore + 
                         Oil_Nonfederal_Offshore + Oil_Federal_Offshore,
                       Gas = Gas_Nonfederal_Onshore +  Gas_Federal_Onshore + 
                         Gas_Nonfederal_Offshore + Gas_Federal_Offshore)]

bh = fread(data.dir%&%'/Baker Hughes/Baker_Hughes_OG_split.csv')
setnames(bh, tolower(names(bh)))
bh[, date:=as.Date(date)]


max.date = as.Date('2020-06-01')

pdf(output.dir%&%'/Model_validation_'%&%today()%&%'.pdf', family='serif', width=11, height=8.5)
par(mai=par('mai')*c(1,1,1,2))
# par(mai=par('mai')*c(1,1,1,2.5))
plot(Oil ~ date, data=spud.sim[date<=max.date], type='n', #lwd=3, col='#74645e',
     xlim = c(spud.sim[, min(date)], as.Date('2020-06-30')),
     ylim=c(0,3000), ylab='Wells drilled per month', xlab='', cex.axis=1.5, cex.lab=1.5)
# Actuals
lines(Oil ~ date, data=spud.sim[date<=sim.start.date], type='l', lwd=3, col='#74645e')
lines(Gas ~ date, data=spud.sim[date<=sim.start.date], type='l', lwd=3, col='#88c4f4')
# Simulated
lines(Oil ~ date, data=spud.sim[between(date, sim.start.date, max.date)], type='l', lwd=3, col='#74645e', lty=5)   
lines(Gas ~ date, data=spud.sim[between(date, sim.start.date, max.date)], type='l', lwd=3, col='#88c4f4', lty=5)
abline(v=seq.Date(as.Date('2010-01-01'), as.Date('2021-01-01'), by='2 years'),
       lty=3, col='lightgray')
abline(v=sim.start.date, lty=1, col='gray50', lwd=2)
text(x=sim.start.date, y=2900, labels='Historical\nPeriod', pos=2, cex=1.5)
text(x=sim.start.date, y=2900, labels='Simulation\nPeriod', pos=4, cex=1.5)

legend(x=as.Date('2015-03-01'), y=2600, lwd=3, lty=1, col=c('#74645e','black'), bty='n',
       y.intersp = 0.75, cex=1.25,
       legend=c('Oil wells drilled (left)', 'Oil rigs active (right)'))

legend(x=as.Date('2011-06-01'), y=350, lwd=3, lty=1, col=c('#88c4f4','blue'), bty='n',
       y.intersp = 0.75, cex=1.25,
       legend=c('Gas wells drilled (left)', 'Gas rigs active (right)'))

legend(x=as.Date('2018-11-15'), y=1800, lwd=2, lty=2, col='gray20', bty='n', cex=1.2,
       legend=c('Simulated\nwells\ndrilled'))

par(new=TRUE)
plot(oil ~ date, data=bh[between(date, spud.sim[, min(date)], max.date)], type='l', lwd=3,
     xlim = c(spud.sim[, min(date)], as.Date('2020-06-30')),
     col='black', ylim=c(0,2000), yaxt='n', xaxt='n', ylab='', xlab='')
lines(gas ~ date, data=bh[between(date, spud.sim[, min(date)], max.date)], type='l', lwd=3, col='blue')
axis(4, cex.axis=1.5)
mtext('Rigs active', side=4, line=2, cex = 1.5)

dev.off()
